close all
clear all
clc
%%
% % =========================
Eng=8.05;
lambda=1.54;
Attl=1008e4; % Attenuation length in unit of Angstrom for water @ 8.05 keV
mu=1/Attl;
beta=mu*lambda/4/pi;

I0=1e2;  
% % =========================
Fluores_Fe_Qz_file={
    'SMFe01_175to212_Fe_EL.txt'; % Fe & 100mM KCl
    'SMFe01_221to258_Fe_EL.txt'; % Fe & 50mM KCl
    'SMFe01_267to304_Fe_EL.txt'; % Fe & 10mM KCl
    'SMFe01_313to350_Fe_EL.txt'; % Fe & 1 mM KCl
    'SMFe01_44to81_Fe_EL.txt';  % Fe No KCl
    
    'SMFe01_30to35_Fe_EL.txt';  % Pure water subphase
    'SMFe01_451to476_Fe_EL.txt'; % Pure 4 mM FeCl3 
    
    }

Fluores_DB_Qz_file={
    'SMFe01_175to212_DB_EL.txt'; % Fe & 100mM KCl
    'SMFe01_221to258_DB_EL.txt'; % Fe & 50mM KCl
    'SMFe01_267to304_DB_EL.txt'; % Fe & 10mM KCl
    'SMFe01_313to350_DB_EL.txt'; % Fe & 10mM KCl
    'SMFe01_44to81_DB_EL.txt'; % Fe No KCl
    
    'SMFe01_30to35_DB_EL.txt';  % Pure water subphase
    'SMFe01_451to476_DB_EL.txt'; % Pure 4 mM FeCl3 
    
    }

qz=[0.01:0.001:0.030]';

MarkerType={'o','s','v','^','d','>','<'};
COLOR={'k','r','b',[0 0.4 0],'m','k','k'};
%% Check direct beam profiles
close all
figure
hold on
for k=1:length(Fluores_DB_Qz_file)
    data=dlmread(Fluores_DB_Qz_file{k})
    DBEL{k}.Intensity=data(:,2);
    DBEL{k}.Intensity_err=data(:,3);
    H=ploterr(qz,DBEL{k}.Intensity,[],DBEL{k}.Intensity_err,MarkerType{mod(k-1,7)+1},'hhy',0.2)
    set(H,'LineWidth',3,'MarkerSize',8);
end
hold off
%set(gca,'yscale','log');
%%
close all
figure
hold on
DL=length(Fluores_Fe_Qz_file)

for k=1:length(Fluores_Fe_Qz_file)
    data=dlmread(Fluores_Fe_Qz_file{k});
    FeEL{k}.Intensity=data(:,2);
    FeEL{k}.Intensity_err=data(:,3);
    H=ploterr(qz,FeEL{k}.Intensity,[],FeEL{k}.Intensity_err,MarkerType{mod(k-1,7)+1},'hhy',0.2)
    set(H,'LineWidth',3,'MarkerSize',8,'color',COLOR{mod(k-1,7)+1});
end
hold off
set(gca,'LineWidth',3,'box','on');

Intensity_sub=FeEL{DL}.Intensity-FeEL{DL-1}.Intensity;
Intensity_sub_err=sqrt(FeEL{DL}.Intensity_err.^2+FeEL{DL-1}.Intensity_err.^2);


%%
MarkerType={'o','s','v','^','d','>','<'}
COLOR={'k','r','b',[0 0.4 0],'m','k','k'}

% Subphase only profile fitting
close all
clc
Param=[
    1, 0.0219, 0
    ]
LB=[0, 0.0218-1e-6,  0];
UB=[1, 0.022+1e-6,  1 ];

xfit=(0.01:0.0001:0.0300)';
xdata_sub=qz(Intensity_sub>0 & qz<0.031);
ydata_sub=I0*Intensity_sub(Intensity_sub>0 & qz<0.031);
ydata_sub_err=I0*Intensity_sub_err(Intensity_sub>0 & qz<0.031);
[p_opt_sub]=optim_fitsubphase(Param,xdata_sub,ydata_sub,ydata_sub_err,LB,UB)
C_bulk=p_opt_sub(1);
Qc=p_opt_sub(2);
bkg=p_opt_sub(3);
T=Transmission(xfit,Qc,beta,lambda);
D=fitsubphase(xfit,Qc,Eng, Attl);
yfit_sub=C_bulk*T.*D+bkg;
figure
hold on
H=ploterr(xdata_sub,ydata_sub,[],ydata_sub_err,MarkerType{1},'hhy',0.1);
set(H,'color',COLOR{1},'LineWidth',3,'MarkerSize',8);
plot(xfit,yfit_sub,'-','LineWidth',3,'color',COLOR{k});
hold off
set(gca,'Linewidth',3,'box','on');
set(gca,'FontSize',30,'FontName','times');
set(gca,'xlim',[0.009,0.036]);
% set(gca,'ylim',[-0.1 3.1]);
axis square

%%
close all
clc


figure
hold on
for k=1:5
    FeEL_corr{k}.Intensity=I0*(FeEL{k}.Intensity-FeEL{DL-1}.Intensity);

    FeEL_corr{k}.Intensity_err=I0*sqrt(FeEL{k}.Intensity_err.^2+FeEL{DL-1}.Intensity_err.^2);
    
    H=ploterr(qz,FeEL_corr{k}.Intensity,[],FeEL_corr{k}.Intensity_err,MarkerType{mod(k-1,7)+1},'hhy',0.1);
    set(H,'LineWidth',3,'color',COLOR{mod(k-1,7)+1})
end
hold off

%%
close all
clc

figure
hold on
xfit=(0.01:0.0001:0.0300)';
qmin=1;
qmax=21;
C_EL=zeros(1,5);
Zion=zeros(1,5);
for k=1:5
    x=qz;
    y=FeEL_corr{k}.Intensity;
    yerr=FeEL_corr{k}.Intensity_err;
    
    xdata=x(qmin:qmax);
    ydata=y(qmin:qmax);
    ydata_err=yerr(qmin:qmax);
    
    
    %%%%%%%%%%%%%%%%%%%%% Radiation Damage Model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Param=[
%         1, 0.0218, 0, 0.1,0,0.021
%         ]
%     LB=[0, 0.0218,  20, 0, 0.022];
%     UB=[100, 0.0220, 100,1,0.035];
%     [p_opt{k}]=optim_fitsurface_tau(Param,xdata,ydata,ydata_err,LB,UB) 
%     C_EL=p_opt{k}(1);
%     Zion=p_opt{k}(3);
%     tau=p_opt{k}(4);
%     qtau=p_opt{k}(5);
%     
%     T=Transmission(xfit,p_opt{k}(2),beta,lambda);
%     D=fitsubphase(xfit,p_opt{k}(2),Eng, Attl);    
%     yfit{k}=C_EL*T.*exp(-Zion./D).*exp(-0.5*(1+erf((xfit-qtau)/0.0001))*tau.*(xfit-qtau)./qtau);
    
    %%%%%%%%%%%%%%%%%%% Regular Model for esitimation maximum binding %%%%%%%%%%%%%%%%%%%%
    Param=[
        1, 0.0218, 10
        ]
    LB=[0, 0.0217-1e-6,  10];
    UB=[100, 0.022+1e-6, 50];
    
    p_opt=optim_fitsurface(Param,xdata,ydata,ydata_err,LB,UB);
    C_EL(k)=p_opt(1);
    Qc=p_opt(2);
    Zion(k)=p_opt(3);
    T=Transmission(xfit,p_opt(2),beta,lambda);
    D=fitsubphase(xfit,p_opt(2),Eng, Attl);    
    yfit{k}=C_EL(k)*T.*exp(-Zion(k)./D);
    plot(xfit,yfit{k},'-','LineWidth',3,'color',COLOR{k});
    H=ploterr(x,y,[],yerr,MarkerType{k},'hhy',0.1);
    set(H,'color',COLOR{k},'LineWidth',3,'MarkerSize',12,'MarkerFaceColor','w');
  
end
set(gca,'Linewidth',3,'box','on');
set(gca,'FontSize',30)%,'FontName','times');
set(gca,'xlim',[0.009,0.031]);
set(gca,'ylim',[-0.1 35.1]);
axis square


%% Calculate the surface density
% subphase concentratin is 5 mM equivalent to 5*
nb=4*1e-3*6.02e23/1e27;
NumFeDensity=zeros(5,1);
for k=1:5
    NumFeDensity(k)=C_EL(k)/C_bulk*nb;
end
NumFeDensity
Zion
concentration=[2, 5, 10, 20, 40];
close all
clc
figure
plot(concentration, NumFeDensity,'ko','MarkerSize',10,'LineWidth',3);
set(gca,'linewidth',3);
set(gca,'box','on');
set(gca,'FontSize',30);
set(gca,'xlim',[-1 50]);
axis square

% # of Molecule per surface area =  C_surf/C_bulk*nb
%% Display
close all
figure
hold on
plot(xfit,yfit_sub,'-','LineWidth',3,'color','k');
H=ploterr(xdata_sub,ydata_sub,[],ydata_sub_err,MarkerType{1},'hhy',0.1);
set(H,'color','k','LineWidth',3,'MarkerSize',14,'MarkerFaceColor','w');
plot(xfit,yfit{1},'-','LineWidth',3,'color','r');
H=ploterr(xdata,ydata,[],ydata_err,MarkerType{2},'hhy',0.2);
set(H,'color','r','LineWidth',3,'MarkerSize',14,'MarkerFaceColor','w');
hold off
set(gca,'xlim',[0.0095 0.0305],'ylim',[-0.5 35.5]);
set(gca,'fontsize',30,'linewidth',3);
set(gca,'box','on');
axis square


%% Error Analysis
% Comment: due to the potential radiation damage issue, the depth of the
% ion layer is considered to be the main error/uncertainty source
% So just shake up the Zion
close all
clc

Zion_test=[25:1:55];

UB_test=UB;
LB_test=LB;
Param=[
    1, 0.0218, 10
    ]
LB=[0, 0.0219-1e-6,  0];
UB=[100, 0.0219+1e-6, 50];

deltaZion=zeros(1,5);
ZionOpt=zeros(1,5);
chisq_min=zeros(1,5);
IronSurfDensity_opt=zeros(1,5);
IronSurfDensity_err=zeros(1,5);
 
for filenum=1:5
    x=qz;
    y=FeEL_corr{filenum}.Intensity;
    yerr=FeEL_corr{filenum}.Intensity_err;
    
    xdata=x(qmin:qmax);
    ydata=y(qmin:qmax);
    ydata_err=yerr(qmin:qmax);
    
    Z=[];
    chisq=[];
    IronSurfDensity=[];
    for k=1:length(Zion_test)
        UB_test(3)=Zion_test(k)+1e-6;
        LB_test(3)=Zion_test(k)-1e-6;
        [p_test,fval]=optim_fitsurface(Param,xdata,ydata,ydata_err,LB_test,UB_test);
        
        Z=[Z,Zion_test(k)];
        chisq=[chisq,fval];
        IronSurfDensity=[IronSurfDensity, p_test(1)/C_bulk*nb];
        
        T=Transmission(xfit,p_test(2),beta,lambda);
        D=fitsubphase(xfit,p_test(2),Eng, Attl);
        yfit{k}=p_test(1)*T.*exp(-Zion_test(k)./D);
      
    end
   
    figure
    subplot(1,3,1);
    hold on
    coef=polyfit(Z,chisq,2);
    chisq_min=coef(3)-coef(2)^2/coef(1)/4;
    chisqfit=polyval(coef,Z);
    plot(Z,chisq,'ko','MarkerSize',10);
    plot(Z,chisqfit,'k-','linewidth',3);
    hold off
    deltaZion(filenum)=sqrt(1/coef(1))*sqrt(chisq_min) % Rescale the uncertainty in parameter.
    ZionOpt(filenum)=-coef(2)/coef(1)/2
    set(gca,'box','on','LineWidth',3);
    grid on
    subplot(1,3,2);
    hold on
    plot(Z,IronSurfDensity,'bs','MarkerSize',10);
    %IronSurfaceDensity_opt(filenum)=IronSurfDensity( ZionOpt(filenum)-deltaZion(filenum)<Z<ZionOpt(filenum)+deltaZion(filenum))
    coef=polyfit(Z,IronSurfDensity,1);
    IronSurfDensityFit=polyval(coef,Z);
    plot(Z,IronSurfDensityFit,'k-','linewidth',3);
    IronSurfDensity_opt(filenum)=polyval(coef,ZionOpt(filenum))
    IronSurfDensity_err(filenum)=coef(1)*deltaZion(filenum)
    hold off
    set(gca,'box','on','LineWidth',3);
    grid on
   
    subplot(1,3,3);
    hold on
    H=ploterr(xdata,ydata,[],ydata_err,'ko');
    set(H,'MarkerSize',12,'LineWidth',3);
    for k=1:length(Zion_test)
        plot(xfit,yfit{k},'k-','LineWidth',3);    
    end
    hold off
    set(gca,'xlim',[0.009,0.036]);
    set(gca,'box','on','LineWidth',3);
    grid on
end
%%
IronSurfDensity_opt
IronSurfDensity_err
c_KCl=[100,50,10,1,1e-5]
close all
figure
H=ploterr(c_KCl,IronSurfDensity_opt,[],IronSurfDensity_err,'ko','hhy',0.2)
set(H, 'LineWidth',3,'MarkerSize',14)
%set(gca,'xscale','log');
set(gca,'ylim',[-0.002,0.052],'ytick',[0:0.01:0.05]);
set(gca,'xlim',[-5,120],'xtick',[0:20:100]);
set(gca,'linewidth',3,'box','on','FontSize',30);
axis square




% %% plot error analysis results
% close all
% clc
% figure
% subplot(1,2,1)
% %hold on
% H=ploterr(concentration, IronSurfDensity_opt,[],IronSurfDensity_err,'ko','hhy',0.1)
% set(H,'MarkerSize',10,'LineWidth',3);
% set(gca,'linewidth',3);
% set(gca,'box','on');
% set(gca,'FontSize',30);
% set(gca,'xlim',[-1 50]);
% axis square
% 

% % %%
% % 
% % haxes1=gca;
% % 
% % haxes1_pos = get(haxes1,'Position'); % store position of first axes
% % haxes2 = axes('Position',haxes1_pos,...
% %               'XAxisLocation','top',...
% %               'YAxisLocation','right',...
% %               'Color','none');
% %           
% % H=ploterr(concentration,ZionOpt,[],deltaZion,'bs','hhy',0.1)
% % set(H,'MarkerSize',10,'LineWidth',3,'parent','haxes2');
% % % set(gca,'linewidth',3);
% % % set(gca,'box','on');
% % % set(gca,'FontSize',30);
% % % set(gca,'xlim',[-1 50]);
% % % set(gca,'ylim',[-1 51]);
% % % set(H,'yaxislocation','right');
% % axis square
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
